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r  A  mixed-type  (hyperbolic,  parabolic,  elliptic)  partial  differential  equation,  which  results  from  a 
model  of  a  towed  cylindrical  acoustic  antenna  in  the  ocean,  has  been  solved  by  combining  the  techniques  Q 
of  the  method  of  lines  and  the  numerical  ordinary  differential  equation  method.  The  model  describes  the  . 
transverse  dynamic  motion  of  a  zero  bending  rigidity  cylinder  in  an  axial  viscous  stream.  This  initially  _ 

ill-posed  problem  was  satisfactorily  solved  bv  the  above  approach.  The  efficiency  of  this  numerical  solution 

has  been  improved  over  previous  solutions.  Discussions  of  these  improvements  are  the  main  theme  of  this  _ 

paper  which  presents  an  explicit  finite-difference  scheme;  the  improvements  increase  computation  speed  f 
while  maintaining  the  required  accuracy.  The  improvements  are  discussed  both  theoretically  and  — — 

computationally  with  respect  to.  consistency,  stability  and  ^convergence.  The  results  of  a  numerical  test  f  CodSI 


example  are  given^  y\s  1  S 
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1.  INTRODUCTION 
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A  partial  differential  equation  (PDE)  describing  the  transverse  dynamics  of  a  long  cylinder  in  axia 


T 


flow  was  first  obtained  by  Paidoussis  [1],  This  dynamic  motion  equation  is  a  fourth-order  PDE. 
In  many  problems  of  interest  the  aspect  (length-to-diameter)  ratio  is  quite  large  and  the  excitation 
forces  are  much  larger  than  the  diameter.  In  these  cases  the  flexural  stiffness  of  the  cable  is  negligible 
relative  to  the  tension  stiffness  term.  Neglecting  this  term  results  in  a  mainstream  equation  which 
is  a  mixed  type;  varying  from  a  hyperbolic  equation  to  a  parabolic  equation  elliptic  equation.  Over 
a  decade  ago,  a  solution  to  this  problem  was  proposed  by  Ortloff  and  Ives  [2]  using  the  method 
of  separation  of  variables.  Their  solution  contained  Bessel  functions  whose  order  and  argument 
are  complex  and  of  a  magnitude  not  well  approximated  by  asymptotic  solutions.  Moreover,  the 
treatment  of  the  free-end  boundary  condition  was  through  a  “bounded"  criteria  which  is  not 
applicable  to  numerical  solutions.  The  Ortloff  and  Ives  solution  [2]  is  also  quite  inefficient  for 
arbitrary  excitation  time  services  because  the  solution  must  be  repeated  for  each  frequency  present 
in  the  Fourier  services  of  the  excitation  term.  Recently,  Lee  and  Kennedy  [3]  made  the  problem 
numerically  well-posed  by  introducing  an  appropriate  treatment  to  the  free-end  boundary  condition. 
Then,  a  combination  of  the  method  of  lines  and  a  numerical  ordinary  differential  equation  (ODE) 
method  (generalized  Adams  methods)  were  applied  to  solve  the  well-posed  problem.  Details  of  the 
theoretical  and  numerical  developments  can  be  found  in  Ref.  [3],  The  efficiency  of  the  above 
solution  can  be  further  improved  by  a  finite-difference  scheme.  The  main  theme  of  this  paper  is  an 
improved  solution  using  a  basic  explicit  finite-difference  scheme.  Further,  this  development  is 
proved  theoretically  convergent.  Prior  to  the  convergence  proof,  the  theory  with  regard  to  stability 
is  fully  discussed.  The  consistency  condition  is  examined.  The  validity  of  the  finite-difference  solution 
is  demonstrated  by  an  exact  solution  test  which  is  used  to  examine  the  accuracy  of  the  mathematical 
formulation  as  well  as  computational  procedures. 


2.  PROBLEM  AND  SOLUTION  SUMMARY 

A  fourth-order  PDE  fl]  representing  the  dynamic  motion  of  a  cylinder  in  a  viscous  axial  flow 
form  is  as  follows: 
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d 

dtcx 

dx 

t( 

\U2{L-x)C-± 

2 

cx_ 

EI^  +  (M  +  +  MV2pj  +  2  MU 

cx *  dr  dx2 


The  parameters  appearing  in  equation  (1)  have  the  following  definitions: 

El  =  bending  rigidity. 

M  =  the  effective  mass  of  fluid  ''pushed’*  by  the  cylinder  per  unit  length  of  cylinder, 
m  =  mass  of  the  cylinder  per  unit  length, 

U  =  fluid  velocity, 

cT  =  longitudinal  drag  coefficient  due  to  hydrodynamic  forces  acting  along  the  cylinder, 
L  =  total  length  of  the  cylinder 

and 

cN  =  normal  drag  coefficient  due  to  hydrodynamic  forces  acting  normal  to  the  cylinder. 
After  neglecting  the  flexural  rigidity  and  rearranging  terms,  equation  (1)  is  reduced  to 


(1) 


c2y  M  +  m  d2y 
dx2  M  cx2 


U 2 


1  V 


I 'J 2  cx  U2  c2x  1  cx 

mIL-"\  +  £w(ci  +  c*)+2vd-,  +  2c*T;-0  l2) 


Ortloff  and  Ives  [2]  used  the  non-dimensional  terms 

t  =  U/L)U ,  P  =  M/(M  +  m).  <f  =  x/L,  £  =  L/D  and  rj  =  y/L. 
Substitution  of  these  terms  into  equation  (2)  yields 


0  + 


,  (ft  +  *N ,  1  „  n 

rr  +  ~ 2~ £"«  +  2  c"£^  =  ° 


(3a) 


The  dynamic  motion  equation  problem  is  described  by  equation  (3a).  The  associated  initial  and 
boundary  conditions  are: 


r]{ 0,  £)  =  q,(<f),  describing  initial  deflection: 
=  0,  zero  initial  velocity. 


dx 


i  =  0 


t,0)  =  0,  fixed-end  condition; 

1)|  <  oc,  bounded  free-end  deflection. 


(3b) 

(3c) 

(3d) 

(3e) 


The  set  of  equations  (3a— e)  poses  an  initial-boundary  problem  in  the  region  of  interest  (Fig.  1). 

Using  the  method  of  characteristics  we  may  classify  the  solution  types  by  evaluating  the  char¬ 
acteristic  type.  To  do  this  we  evaluate  the  quadratic  equation 


>0  s  <  s*  hyperbolic 
cTe(l  -  s)  =0  c  =  <f*  parabolic 
<0  if  >  <f*  elliptic. 


for  some  0  ^  if*  <  I.  The  point  to  be  noted  is  that  over  the  length  of  the  physical  cylinders  all 
three  characteristic  types  are  found.  The  above  summarizes  the  problem  background,  the  region 
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f 

Fig.  1.  Region  of  interest. 


of  interest  and  the  change  of  equation  type. 

Two  solutions  exist  for  equation  (3a).  One  is  the  Bessel  function  solution,  obtained  by  means  of 
the  separation  of  variables  [2].  This  solution  is  inefficient  for  arbitrary  temporal  excitation  of  the 
cylinder  structure.  Lee  and  Kennedy  [3]  introduced  the  numerical  solution  to  remove  the  above 
inefficiency.  Lee  and  Kennedy  introduced  the  change  of  variables 

t  =  r  (4) 

and 

x  =  x(£,t)  =  0t  +  (1  -4),  (5) 

which  brings  equation  (3a)  into  an  equivalent  form: 

4  +  d(x,  tftjxx  +  Si\x  +  crj,  =  0,  (6a) 

where 

d(x,r)  =  /?J^1  -  ^cTe(x  -  )?r)J  -  /?2, 

B(x,t)  =  ~  efi[cNp  ~  (cT  +  cN)], 
c  =  ^CN£0 

and 

r)(x,t)  =  rj(Pz  +  (1  -  <J),t). 

The  associated  initial  and  boundary  conditions  are: 


r)(x,0)  =  7,(1  -  x),  prescribed  initial  deflection;  (6b) 

>}i(x.0l«-o  =  Pfahl*  ~  *)»  initial  velocity;  (6c) 

r)(x,t)  =  f(t)  at  x  =  (it  +  1,  prescribed  fixed-end  condition;  (6d) 

| r)(x,  t )|  <  x  at  x  =  /?r,  bounded  free-end  condition.  (6e) 


The  set  of  equations  (6a-e)  defines  an  initial-boundary  problem  counterpart  to  equation  (3a). 
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The  ODE  solution  to  equation  (6a),  proposed  by  Lee  and  Kennedy  [3],  uses  the  method  of 
lines  to  discretize  partial  derivatives  with  respect  tc  x  by  central  and  forward  differences.  This 
discretization  formulates  the  PDE,  equation  (6a).  into  a  system  of  second-order  ODEs  in  the 
variable  r,  i.e. 


(*?«)«  +  <*('?, )„  + 


■jcfi  +  S 


+ 


a(*m-i,f)  ,  cTe/?  +  S~\  . 

P  +  h  J,”_1 


(7) 


where  h  =  Ax. 

Then,  another  change  of  variable,  =  H-'m,  brings  equation  (7)  into  a  system  of  first-order 
ODEs  in  the  form 


f\  -  A(x,t)>j  +  g(x,f,»7),  (8) 

where  the  vector  g  contains  two  non-zero  components  which  supply  the  fixed-end  and  free-end 
boundary  conditions. 

Because  the  matrix  A(.x.  r)  turns  out  to  be  independent  of  f,  the  numerical  solution  to  equation 
(8)  with  associated  conditions  can  be  obtained  explicitly  by  means  of  the  generalized  Adams  method 

[3]: 


e?„+1  =  eA',f  +  fi(A/irl(eu  -  0)gn,  (9) 

which  summarizes  the  numerical  solution  to  problem  (3a). 

Note  that  the  computation  of  equation  (9)  is  straightforward  provided  that  the  matrix  exponential 
eA*  can  be  carried  out  simply.  An  accurate  and  economical  computation  of  eA'1  is  by  means  of  a 
second-order  diagonal  Pade  approximation: 


e-"l  u  +  1  A/i  +  ^  (Aft)* 


.am 


3  -  ^  Aft  +  — (A/i): 


(10) 


Even  though  the  eA*  need  be  calculated  only  once,  it  is  necessary  at  every  time  step  to  multiply 
the  matrix  by  a  vector.  For  a  matrix  of  order  M,  it  needs  at  least  2 M 3  operations  initially  to  obtain 
the  approximation  for  the  eA\  then  at  every  time  step' there  is  another  M2  operations  needed  to 
produce  the  result.  In  view  of  the  cost  of  the  computation,  in  both  speed  and  computer  storage, 
we  propose  an  explicit  finite-difference  solution  to  problem  (6a).  We  discuss  the  complete  finite- 
difference  development  in  the  following  sections. 


3.  A  STABLE  FINITE-DIE FERENCE  MODEL 

In  this  section  we  describe  the  complete  development  of  a  finite-difference  scheme  for  the  solution 
of  problem  (6a).  It  begins  with  the  formulation.  Following  the  formulation,  the  theory  with  regard 
to  consistency,  stability  and  convergence  is  presented  for  the  proposed  finite-difference  scheme. 

3.1.  An  Explicit  Finite-difference  Formulation 
To  maintain  generality,  let  us  consider  the  non-homogeneous  equation 


q„  +  dfjxx  -f  fit),  +  cf),  =  Fix.  M. 


(ID 


An  ocean  dynamic  motion  mode! 


79 


where  F(x,t)  physically  models  a  FORCING  function.  Pursuing  the  finite-difference  discretization 
of  equation  (11),  let  us  derive  a  partial  differentiation  expression  for  dfjxz: 

=  (M  +  drjx)x 

=  dxxrj  +  2  dxr)x  +  drjxx.  (12) 


=  /?j^l  -  ^cTe(l  +  x  -  0r)J  -  p2. 


therefore  equation  (12)  becomes 


which  implies  that 


=  ~  2  ct£0 


<*„  =  0, 


(dt?)„  =  drjxx  -  cTeptjx, 


dijxx  =  (<3  r))xx  +  cT£prjx. 


Now,  apply  the  conventional  finite-difference  discretization  to  equation  (11),  making  use  of  equation 
(13),  using  the  central  difference  for  the  second  derivative  and  the  backward  difference  for  the  first 
derivative,  we  obtain  a  different  equation  for  equation  (11),  viz. 

f\v  1  -  2r)”m  +  fc-  1  d”m+  M+i  -  Wm  + 
k2  h2 

An  _  An  An  _  An  ja»  +  1  _  -an 

-  cTef} +  6nm  =  F(xm, t"),  (14) 

h  n  k 

where  k  =  At,  h  =  Ax. 

A  rearrangement  of  equation  (14)  will  result  in  the  simple  expression 


»r 1  +  prm  +  qc-  i  +  +  rm~l  =  sF(xm,n 


where 


AY  k2  ~k2  1 

\h) 

-  k(J\  +  k£). 

(16) 

[1  +M), 

(17) 

r/i,v  i 

=  LW  +  fc<?) 

(18) 
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S  =  k2/(  1  +  kc). 


pv 1  -  sF(xm.n  -  /*•?:  -  qpi-  .  -  rpi.  ,  -  (20a) 

Equation  (20a)  is  the  explicit  finite-difference  scheme  for  the  solution  of  equation  (11).  Initially,  to 
compute  the  the  associated  initial  and  boundary  conditions  (6b— e)  must  be  fulfilled : 

pnm~l  =  ^1(1  -  *J.  the  prescribed  initial  deflection;  (20b) 

r/”  is  obtained  making  use  of  the  initial  velocity  information,  such  that 

Pm  =  Pm~l  +  biPj-nA  l-*);  (20c) 

Pm-i  —  f(0  at  x  =  fit  +  1,  the  prescribed  fixed-end  condition;  (20d) 

pnm  +  i  <  qo  at  xm  +  ,  =  fit,  the  bounded  free-end  condition.  (20e) 

Problem  (20a)  with  associated  conditions  (20b-e)  constitutes  a  counterpart  initial-boundary  problem 
to  problem  (6a). 

3.2.  Consistency 

In  examining  whether  the  difference  equation  (15)  is  consistent  with  the  differential  equation  (11), 
we  expand  p”m* 1  and  p„+l  in  powers  of  k  upon  p"m,  and  expand  pmm+l,pnm.  1  in  powers  of  h  upon 

PI¬ 

VOT  rotational  economy  we  drop  the  indices  n  ?rd  m  so  that  p  =  p1^.  Then 


Pm*  1  =  P  +  Wi  +  ^P"  +  "'■> 


Pm  *  =  P  ~  kp,  +  p„ 


P"m*i  =  P  +  PPx  +2|  P**  + 


Pm- l=p-PPx  +  2jPx*  +  -  - 


To  examine  the  consistency,  we  localize  the  coefficients  <3j  in  the  closed  interval  [t,  t  +  At]  using  a 
to  represent  the  average  We  rearrange  equation  (14)  and  put  all  terms  on  the  l.h.s.  of  the 
equation  to  give 

h2(PV '  -  2P'm  +  P”m  ')  +  k2a{P'm.  ,  -  2P"m  +  PI. ,) 

-  (crefi  -  bWm  ~  Pm- .)*:*  +  ckh2iPV  1  -  PI)  -  k2h2F[xm.t )  =  0. 


(25) 
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Substituting  equations  (2 1 )— (24)  into  equation  (25),  we  obtain 
k2h2fj„  +  k2h2arjxx  -  (cTe/?  -  i>)k2h2fjx 

+  jj^rj,,,,  +  k2h*arjxxxx  +  (cTe/?  -  b)k2h3rjx 

+  ck2h2r)y  +  ^ k3h2fj„  -  k2h2F(xm,t)  =  0. 


k2h2[n,<  +  (arjxx  -  crepr)x)  +  bf\x  +  crj,  -  F(xm,r)] 


c  *  lc^  k2h* 

=  ~-k3h2rj „  -  (crep  -  b)k2h3flxx  +  —  h2fj,„  -  — afjxxx ; 


the  truncation  error  E  is  given  by 


c  »  k2  h 2 

E  =  ~  j  ~  -  b)hrjxx  +  —»/<«.<  -  -^afixxxx  +  high-order  terms 

=  0(k,h). 


In  the  ^lim^  £  =  0,  which  establishes  the  consistency. 


3.3.  Stability 


From  equation  (14),  if  one  compares  the  terms,  in  equation  (14),  between  ar|d  il 

is  easily  seen  that  as  k  -*0,  the  dominant  term  is  Similarly,  the  dominant  term  between 

o(j^J  a°d  °(^)  ^>e,on8s  to  the  term.  This  analysis  suggests  that  for  the  stability  analysis. 


one  can  use 

fj<,  +  =  0  (27) 

to  derive  the  sufficient  condition  of  stability. 

Theorem 

The  finite-difference  scheme  (14)  used  to  solve  the  dynamic  motion  equation  (1 1)  is  stable  if 

*<± 

h  "  |df 

Proof.  For  brevity,  write  /.  =  e*‘.  Since  d(x)  is  not  rapidly  varying  in  the  variable  x,  we  adopt 
the  “frozen  variable”  approximation  to  localize  <3(x)  as  a  constant  in  x  between  the  interval 
[t,  t  +  At]  and  also  for  brevity,  write  d(x)  =  d  =  const.  We  look  for  the  solution  i)(x,  t)  in  the  form 
e"e‘w*,  using  the  Von  Neumann  approach.  Substituting  =  e*'e,“i  into  equation  (27),  using  the 
approximation  a  =  d(x),  one  finds 


♦  i » _ 2  c*i*  jf.  g**"  -  1  ** 


*  1th  _  2  -  1  >h 


•  e'“"  +  a  ■ 


■  e*"*  =  0. 
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A  simplification  of  the  above  equation  gives 

;.2  +  K?  +  1=0, 


(28) 


where 


K  =  2 


k2 

a— i  [cos(wh)  -  1]  -  1 

n 


Then 


/.  = 


-  K  ±  JK2  -  4 
2 


(29) 


The  stability  condition  requires  that  |/j  S;  1,  which  implies  that  both  of  the  following  conditions 
must  be  satisfied: 


and 


-K  +  jK2-  4 


$  1 


/.3  = 


K-jK2-  4 


^  1. 


(30) 


(31) 


Both  conditions  (30)  and  (31)  have  the  term  S/K2  —  4  in  common.  This  is  the  term  that  needs  to 
be  analyzed  in  order  to  derive  the  stability  condition.  Condition  (30)  is  equivalent  to 


K 

■> 


1 

?V 


(32) 


There  are  two  cases  to  examine: 


(i)  K2  -  4  >  0 


(33) 


and 


(ii>  K2  -  4  <  0, 


(34) 


Case  (i): 


K2  —  A  >  0  -*  K2  >4 


2<jfl^[cos(cu(i)  -  1] 


>  4. 


In  order  that  the  above  condition  is  satisfied,  it  is  necessary  that  [cosfcuh)  -  1]  #  0.  Since  a  <  0 
in  the  hyperbolic  region  and  [cos(cofc)  -  1]  <  0.  together  they  require  that 
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which  implies  that 


u  >  {a[cos(toh)  -  l]}'1. 


Returning  to  condition  (32),  we  examine 


-f +  ^v/K2  -  4<  1  -8, 


which  always  exists. 
Also, 


-I  <  ~f  +  \y/K2 


K  1  r~T 


since  condition  (36)  is  always  satisfied;  which  also  implies  that  -  —  —  -%/ K 2  —  4  ^  1  is  obviously 

true. 

Next, 

-1  <  -  y  ~  |  v/^2  -  4  —  /C  ^  2.  (38) 


On  close  examination  of  conditions  (35)— (38),  it  is  clearly  seen  that  the  case  K2  -  4  >  0  will  not 
lead  to  a  condition  of  stability. 


Case  (ii): 


K2-4<0^  “f  +  \VK2~*  <  1 


-  y  +  ^  V/K2  -  4  =  1  identically. 


Therefore,  K2  ^  4  gives 


a jjj [cos (a»/t)  -  1]  ^  <  2(a)"  1[cos(ct»/*)  —  1]  *, 


max  [cos(ajh)  -  1]  =  2  gives  [equation  (39)] 


k2  ^  1  k  ^  1 

n  ^  r  ^  ~~r' 

h2  a  h 


establishing  the  condition  of  stability. 


Q.E.D. 


3.4.  Convergence 

The  consistency  condition,  established  in  Section  3.2,  states  that  the  finite-difference  solution 
approaches  the  true  solutions  as  Ic,  h  ->  0.  Additionally,  the  stability  condition,  proved  in  Section 
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3.3,  states  that  the  numerical  solution  approaches  the  finite-difference  solution  because  at  every 
time  step  the  numerical  error  remains  bounded.  Both  the  stability  and  consistency  imply  that  the 
numerical  solution  approaches  the  true  solution  as  n  -*  x  and  /c,  /i  — ►  0;  this  is  the  concept  of 
convergence. 


4.  NUMERICAL  COMPUTATIONS 

This  section  presents  an  exact  solution  test  [3]  which  is  used  to  examine  the  theoretical  and 
computational  accuracy  of  the  finite-difference  model.  All  computations  were  performed  on  the 
VAX  11/780  computer  using  double-precision  complex  arithmetic. 

Consider  the  function 


iifi;,  r)  =  e'^sinfcu^).  (41) 

If  we  perform  the  partial  derivations  i;,t  and  ijr,  and  substitute  them  into  the  original 

dynamic  motion  equation  (3a),  we  ohtain  an  inhomogeneous  equation  of  the  form 

f!,r  +■  d&Ui  +  +  2/?,.t  +  CE/r  =  f(C.T).  (42) 


where 


Fl;,  r)  =  *J  -w2[a(^)  +  1]  +  ia>c}  sin  (ox;)  +  (bio  +  2  Piio2)  cos  (u)<;fe,a\  (43) 

We  regard  F(c,  r)  as  simulating  a  forcing  function  and  select  input  parameter  values  close  to  those 
of  realistic  physical  parameters. 

The  associated  initial  and  boundary  conditions  are  taken  as  follows: 

e‘"'sin(cu^)|r  =  0  =  sin(w^),  initial  condition;  (44) 

e'“'sin(cuc)|4  =  0  =  0.  fixed-end  condition;  (45) 

ellu,sin(ax;)|»  =  s.B  =  e’"' sin  (a><|;B),  such  that  je‘w,sin  (a»cB)l  $  1  <  x,  free-end  condition.  (46) 
We  purposely  select  u>  =  i  so  that  as  t  -»  x,  rj(^,T)  -*  0.  In  selecting  a  practical  e  =  5E  +  4,  we  find: 


Problem 

Hyperbolic 

Parabolic 

Elliptic 


Region 

0  sS  c  <  0.998 
S  =  0.998 
0.998  <4  $  1 


Note  that  the  elliptic  region  is  extremely  small  and  considered  negligible.  This  is  usually  true  in 
many  problems  of  interest. 

The  following  inputs  are  selected  for  the  computation: 


0=1/2,  cT  =  0.01,  cN  =  0.09,  r?  =  5E  +  4.  ru  =  0.1i.  At  =  IE -5. 


This  problem  is  chosen  such  that  the  real  part  is  zero  and  calculated  up  to  10,000  steps.  The  results 
are  presented  in  Table  1,  and  show  very  good  agreement  with  the  exact  solution. 

If  these  computations  had  been  carried  out  using  the  numerical  ODE  solution,  they  would  have 
required  calculation  of  a  matrix  exponential,  a  matrix  inversion,  two  matrix  matrix  multiplications 
and  two  matrix-vector  multiplications.  The  total  operation  would  be  0(/V3).  The  improved  finite- 
difference  scheme  only  requires  the  calculation  of  four  scalar  vector  multiplications,  four  vector 
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Table  I  An  exact  solution  test-' 


> 

T 

03052E 

-  04 

0  3052E 

-  01 

0  1 1 14E 

*  00 

0  1 

0  lOOOE 

-  01 

0.9967E 

-  02 

0  9886E 

-  02 

0  I000E 

-  01 

09969E 

-  02 

09889E 

-  02 

0  2 

0.2000E 

-  01 

0 1994E 

-  01 

0.I978E 

-  01 

0  2000 E 

-  01 

0.I994E 

-  01 

0  1977E 

-  01 

0.3 

0  3000E 

-  01 

0  2992 E 

-  01 

02970E 

-  01 

0  3000 E 

-  01 

0.2991  E 

-  01 

0  2967 E 

-  01 

04 

0  4000E 

-  01 

0.3988E 

-  01 

0.3957E 

-  01 

0.400  IE 

-  01 

0  3988 E 

-  01 

03956E 

-  01 

0  5 

0. 5002  E 

-  01 

04986E 

-  01 

0.4946E 

-  01 

05002E 

-  01 

0.4986E 

-  01 

04946E 

-  01 

06 

06003E 

-  01 

0  5985 E 

-01 

0.5937E 

-  01 

06003E 

-  01 

0  598 5 E 

-01 

0  5937E 

-  01 

0  7 

0  7005 E 

-  01 

0.6984E 

-  01 

0  6928E 

-  01 

0.7005E 

-  01 

0.6984E 

-  01 

0.6928E 

-  01 

0.8 

0.8008E 

-  01 

07984E 

-01 

07920E 

-  01 

08008E 

-  01 

0.7984  E 

-  01 

07919E 

-  01 

09 

0.9012E 

-  01 

08990E 

-  01 

0.8925E 

-  01 

0.9012E 

-  01 

0.8984E 

-  01 

0.89 12E 

-  01 

'The  table  entries  show  both  the  exact  (upper  xaluel  and  the  numerical 
(lower  value!  solution 


vector  additions  and  a  few  scalar-scalar  operations.  The  total  operation  is  0(N2).  It  is  evident  that 
the  improved  finite-difference  solution  is  much  faster. 


5.  CONCLUSIONS 

The  previous  solution  [3]  to  a  well-posed  dynamic  motion  equation  was  made  available  using 
an  efficient  numerical  ODE  solution  to  a  vibration  problem  of  practical  interest.  This  numerical 
solution  improved  the  generality  and  computational  efficiency  of  the  previous  special  function 
solution.  This  efficient  solution  is  now  further  improved  by  a  stable  explicit  finite-difference  scheme 
which  increased  the  computation  speed  while  maintaining  the  same  accuracy.  The  theoretical 
development  assures  the  stability,  consistency  and  convergence  of  the  solution.  One  test  example 
was  calculated  and  the  answers  compared  very  favorably  with  the  exact  solution,  demonstrating 
the  accuracy.  This  improved  model,  thus  far,  appears  to  give  a  more  efficient  numerical  treatment 
of  the  transverse  dynamics  of  a  long  thin  cylinder  in  viscous  axial  fluid  flow. 
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